function [min_male, max_male, min_female, max_female, min_ceb_male, max_ceb_male, min_ceb_female, max_ceb_female, min_intrashare_male, max_intrashare_male, problem] = get_ricebs(couple,singlemale,singlefemale,leisurem,leisuref,hworkm,hworkf,wagem,wagef,q,Q,nonlabor,mutually_acceptable,divorce_costsM,divorce_costsF,divorce_costsMF,hcat,text_solver)

% ================
% define constants
% ================

T = 112;                                              % total potential working hours
H = length(leisurem);                                 % total number of households in the data
alpha = 0.4;                                          % bound for non-labor income

% =========================
% define decision variables
% =========================
yalmip('clear')
qM = sdpvar(H,1);                                     % private consumption of Hicksian good male
qF = sdpvar(H,1);                                     % private consumption of Hicksian good female

pM = sdpvar(H,H,'full');                              % personalized public good price male
pF = sdpvar(H,H,'full');                              % personalized public good price female

pmhworkM = sdpvar(H,H,'full');                        % personalized male's hwork price male
pmhworkF = sdpvar(H,H,'full');                        % personalized male's hwork price female

pfhworkM = sdpvar(H,H,'full');                        % personalized female's hwork price male
pfhworkF = sdpvar(H,H,'full');                        % personalized female's hwork price female

wageM = sdpvar(H,1);                                  % shadow wage male
wageF = sdpvar(H,1);                                  % shadow wage female

nl_incomeM = sdpvar(H,1);                             % non-labor income male
nl_incomeF = sdpvar(H,1);                             % non-labor income female

cebm = sdpvar(H,1);                                  % consumption based CEB male
cebf = sdpvar(H,1);                                  % consumption based CEB female

sharem = sdpvar(H,1);                                % consumption based RICEB male
sharef = sdpvar(H,1);                                % consumption based RICEB female

intrasharem = sdpvar(H,1);                           % intrahousehold consumption share male


%%
% ==================
% define constraints
% ==================

Constraints = [];

% restriction on individual private consumption
for i=1:H
    Constraints = [Constraints, q(i) == qM(i) + qF(i), qM(i) >=0, qF(i) >= 0];

    if singlefemale(i) == 1
        Constraints = [Constraints,qM(i) == 0];
    elseif singlemale(i) == 1
        Constraints = [Constraints,qF(i) == 0];
    end

end


% shadow price should be non-negative and equal to observed wage if the individual is employed
for i=1:H

    Constraints = [Constraints,wageF(i) >= 0, wageM(i) >= 0];

    if ~isnan(wagem(i))
        Constraints = [Constraints, wageM(i) == wagem(i)];
    end
    if ~isnan(wagef(i))
        Constraints = [Constraints, wageF(i) == wagef(i)];
    end

end


% restriction on public prices
for i=1:H
    for j=1:H

        Constraints = [Constraints, 1 == pM(i,j) + pF(i,j), pM(i,j) >=0, pF(i,j) >= 0];
        Constraints = [Constraints, wageM(i) == pmhworkM(i,j) + pmhworkF(i,j), pmhworkM(i,j) >=0, pmhworkF(i,j) >= 0];
        Constraints = [Constraints, wageF(j) == pfhworkM(i,j) + pfhworkF(i,j), pfhworkM(i,j) >=0, pfhworkF(i,j) >= 0];

        if (couple(i) == 1 || singlemale(i) == 1) && singlemale(j) == 1
            Constraints = [Constraints, 1 == pM(i,j), wageM(i) == pmhworkM(i,j), wageF(j) == pfhworkM(i,j)];
        end

        if (couple(j) == 1 || singlefemale(j) == 1) && singlefemale(i) == 1
            Constraints = [Constraints, 1 == pF(i,j), wageM(i) == pmhworkF(i,j), wageF(j) == pfhworkF(i,j) ];
        end
    end

end


% feasibility restrictions on the individual nonlabor income (for each matched pair)
% individual nonlabor incomes after divorce must lie between 40% and 60% of total nonlabor income under marriage
for i=1:H

    Constraints = [Constraints, nonlabor(i) == nl_incomeM(i) + nl_incomeF(i)];

    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints, alpha*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= (1-alpha)*nonlabor(i), alpha*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= (1-alpha)*nonlabor(i)];
        else
            Constraints = [Constraints,(1-alpha)*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= alpha*nonlabor(i), (1-alpha)*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= alpha*nonlabor(i)];
        end
    elseif singlefemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
    end

end


% individual rationality
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, wageM(i)*T + nl_incomeM(i) - divorce_costsM(i) <= qM(i) + wageM(i)*leisurem(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
        Constraints = [Constraints, wageF(i)*T + nl_incomeF(i) - divorce_costsF(i) <= qF(i) + wageF(i)*leisuref(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
    end
end


% no blocking pairs
for i = 1:H
    for j = 1:H
        if mutually_acceptable(i,j) == 1
            Constraints = [Constraints, (wageM(i) + wageF(j))*T + nl_incomeM(i) + nl_incomeF(j) - divorce_costsMF(i,j) <= ...
                qM(i) + qF(j) + wageM(i)*leisurem(i) + wageF(j)*leisuref(j) + pM(i,j)*Q(i) + pF(i,j)*Q(j) + ...
                pmhworkM(i,j)*hworkm(i) + pmhworkF(i,j)*hworkm(j) + pfhworkM(i,j)*hworkf(i) + pfhworkF(i,j)*hworkf(j)];
        end

    end
end


% define RICEBs and CEBs
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints,
            sharem(i) == (qM(i) + Q(i))/(q(i) + Q(i)),...
            sharef(i) == (qF(i) + Q(i))/(q(i) + Q(i)),...
            cebm(i) == qM(i) + Q(i),...
            cebf(i) == qF(i) + Q(i),...
            intrasharem(i) == (qM(i) + Q(i))/(q(i) + 2*Q(i))];
    elseif singlemale(i) == 1
        Constraints = [Constraints, sharem(i) == 1, cebm(i) == q(i) + Q(i), intrasharem(i) == 1];
    elseif singlefemale(i) == 1
        Constraints = [Constraints, sharef(i) == 1, cebf(i) == q(i) + Q(i)];
    end
end



%%
% =================
% Solve the problem
% =================

min_male = zeros(16,1);
max_male = zeros(16,1);
min_female = zeros(16,1);
max_female = zeros(16,1);

min_ceb_male = zeros(16,1);
max_ceb_male = zeros(16,1);
min_ceb_female = zeros(16,1);
max_ceb_female = zeros(16,1);

min_intrashare_male = zeros(16,1);
max_intrashare_male = zeros(16,1);

problem = 0;
options = sdpsettings('verbose',0,'solver',text_solver);

% males
j = 1;
for empym = 1:2
    for edum = 1:2
        for empyf = 1:2
            for eduf = 1:2

                k = empym*1000 + edum*100 +  empyf*10 + eduf;
                index = (couple == 1 & hcat == k);

                if sum(index) > 0 && problem == 0

                    obj1 = sum(sharem.*index)/sum(index);
                    obj2 = sum(sharef.*index)/sum(index);
                    obj3 = sum(cebm.*index)/sum(index);
                    obj4 = sum(cebf.*index)/sum(index);
                    obj5 = sum(intrasharem.*index)/sum(index);


                    %--------- RICEB Male --------%
                    sol = optimize(Constraints,obj1,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        min_male(j) = value(obj1);
                    end

                    sol = optimize(Constraints,-obj1,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        max_male(j) = value(obj1);
                    end

                    %--------- RICEB Female --------%
                    sol = optimize(Constraints,obj2,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        min_female(j) = value(obj2);
                    end

                    sol = optimize(Constraints,-obj2,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        max_female(j) = value(obj2);
                    end


                    %--------- CEB Male --------%
                    sol = optimize(Constraints,obj3,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        min_ceb_male(j) = value(obj3);
                    end
                    
                    sol = optimize(Constraints,-obj3,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        max_ceb_male(j) = value(obj3);
                    end

                    %--------- CEB Female --------%
                    sol = optimize(Constraints,obj4,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        min_ceb_female(j) = value(obj4);
                    end

                    sol = optimize(Constraints,-obj4,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        max_ceb_female(j) = value(obj4);
                    end


                    %--------- Intrashare Male --------%
                    sol = optimize(Constraints,obj5,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        min_intrashare_male(j) = value(obj5);
                    end

                    sol = optimize(Constraints,-obj5,options);
                    if sol.problem ~= 0
                        problem = 1;
                    else
                        max_intrashare_male(j) = value(obj5);
                    end

                else
                    min_male(j) = NaN;
                    max_male(j) = NaN;
                    min_female(j) = NaN;
                    max_female(j) = NaN;
                    min_ceb_male(j) = NaN;
                    max_ceb_male(j) = NaN;
                    min_ceb_female(j) = NaN;
                    max_ceb_female(j) = NaN;
                    min_intrashare_male(j) = NaN;
                    max_intrashare_male(j) = NaN;
                end
                j = j + 1;
            end
        end
    end
end





return


